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Section  1 
Introduction 


This  report  describes  the  research  and  development  work  that  was  carried  out  between  August  1, 
1997  and  August  1,  1999  under  the  DSWA  SBIR  Phase  II  Program  in  the  area  of  High-Order 
Godunov  (HOG)  Schemes  for  Multiphase  Gas-Particulate  Flowfield  Modeling  and  Simulation. 

We  describe  here  the  development  and  extension  of  a  major  new  capability  for  the  HOG 
technology  begun  in  our  Phase  I  work  and  continued  here.  The  new  code  developed  here  is  based 
on  both  the  Godunov  scheme  solution  to  the  compressible  Navier-Stokes  equations  and  the 
projection  method  solution  of  the  incompressible  Navier-Stokes  equations  and  provides  high 
accuracy  solutions  to  multi-phase  flows  ranging  from  high  speed  compressible  flows  to  long  time 
scale  buoyant  incompressible  flows.  This  code  has  been  used  to  simulate  the  two-phase  plumes 
resulting  from  both  confined  and  unconfined  explosions  of  conventional  warheads  as  well  as  to 
simulate  the  cloud  rise  phenomena  associated  with  nuclear  bursts  in  dry  climates. 

The  Phase  II  objectives  described  in  our  proposal  can  be  summarized  as  follows:  to  simulate 
gas-particulate  cloud  rise  resulting  from  the  compromise  of  both  confined  and  unconfined  CB 
facilities,  to  analyze  the  results  and  compare  with  experimental  data  and/or  other  models  as 
available,  and  to  modify,  implement,  and  optimize  algorithms,  schemes,  and  codes  in  support  of 
these  simulations. 

We  have  successfully  met  all  of  our  objectives.  A  full  simulation  package,  titled  CATAPULT 
(Computational  All-Time  Algorithm  Package  with  Underlying  Lawrence  Technology)  has  been 
developed  to  model  multi-phase  flows  from  explosion  to  fallout,  spanning  both  compressible  and 
incompressible  fluid  regimes.  Simulations  of  explosions  from  both  confined  and  unconfined  CB 
facilities  have  been  completed  and  the  results  analyzed  and  compared  to  existing  data  and  models 
where  applicable.  Several  modifications  and  optimizations  have  been  included  in  the  CATAPULT 
package  to  provide  faster  and  more  accurate  solutions.  In  addition,  simulations  of  nuclear  cloud 
rise  phenomenon  have  also  been  conducted  using  CATAPULT  at  the  behest  of  the  program 
director. 


The  sections  below  consist  of  a  description  of  the  governing  equations,  physical  models  and  several 
extensions  used  here  (Section  2),  the  numerical  algorithms  and  extensions  of  our  original  codes  by 
Krispin  and  Collins  ([19]  and  Collins  et  al.  [12])  are  given  in  Section  3,  and  Section  4  describes  the 
numerical  results.  The  work  is  summarized  in  Section  5,  including  our  main  conclusions  and 
several  research  and  development  issues  that  we  find  important  and  which  we  intend  to  pursue  in 
future  work.  Included  with  this  report  is  a  compact  disk  containing  the  CATAPULT  source  code 
and  instructions  on  compiling  and  running  the  model.  Numerous  results  and  animations  of  the 
work  performed  here  may  also  be  found  on  the  included  compact  disk  or  at  our  website 
http://www.krispintech.com.  This  work  also  resulted  in  two  publications,  ’’Projection  Methods 
for  Incompressible  Multiphase  Cloud-Rise  Phenomena”  [8]  and  ’’Two  Phase  Plume  Dynamics— a 
Comparison  of  Axisymmetric  and  Three  Dimensional  Simulations”  [24],  The  papers  were  both 
presented  at  the  American  Institute  of  Aeronautics  and  Astronautics  summer  meeting  in  Norfolk, 
VA  and  are  included  in  electronic  form  on  the  enclosed  CD. 
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Section  2 

Governing  Equations 


2.1  Introduction. 

The  work  accomplished  during  this  program  included  modeling  multi-phase  flows  in  both  the 
compressible  and  incompressible  fluid  regimes.  The  code  used  for  the  incompressible  regime  is 
based  upon  the  work  completed  in  Phase  I  [17]  where  underlying  code  technology  used  is  the 
projection  method  for  (single-fluid)  incompressible  flow  that  was  introduced  by  Chorin  [10]  and 
updated  in  the  late  1980’s  by  Bell  et  al.  [5]  by  using  an  unsplit  high-order  Godunov  scheme  to 
discretize  the  convection  terms.  In  past  work,  Krispin  and  Collins  [19]  extended  the  projection 
method  to  two-phase,  dusty-gas  flows  by  incorporating  ideas  from  Collins  et  al.  [12]  into  the 
variable  projection  methodology  of  Bell  and  Marcus  [6].  For  this  work,  the  methodology  of 
Krispin  and  Collins  [19]  was  used  as  the  baseline  code  and  it  was  developed  and  extended  to  the 
point  that  allows  to  simulate  post-explosion,  dusty-gas  cloud  rise  dynamics.  The  compressible 
regime  of  the  flow  is  solved  using  the  technique  of  Collins  et  al.  [12]  and  incorporates  full  adaptive 
mesh  refinement  (AMR)  with  gas  and  dust  phases  coupled  in  the  same  manner  as  in  the 
incompressible  code. 

Reference  [19]  describes  two  numerical  schemes  for  solving  incompressible  gas  flows  with 
embedded  dust  particles.  The  flow  is  modeled  by  extending  the  ideas  described  in  Collins  et  al. 
[12]  for  the  dusty-gas  compressible  fields  to  the  incompressible  gas  limit.  First,  the  Eulerian 
compressible  conservation  laws  are  used  for  the  dust  phase  in  the  dusty-gas  limit,  where  the  dust 
particle  pressure  is  zero  and,  as  a  consequence,  all  the  eigenvalues  are  linearly  degenerate  and  the 
fast  modes  are  absent  (see  for  example  [11],  [22]).  For  the  gas  phase,  the  incompressible  Eulerian 
equations  are  used.  The  two  phases  are  coupled  via  drag  forces  that  are  described  as  sources  in 
the  momentum  equations  of  both  phases.  Also,  gravity  forces  act  on  both  phases.  Following  the 
work  of  [12[,  the  limit  of  the  resulting  equations  for  a  velocity  equilibrium  of  both  phases  was 
analyzed  in  [19];  it  was  found  that  equilibrium  solutions  to  this  system  exist  only  when  the 
pressure  gradient  vanishes,  resulting  in  a  trivial  solution.  The  above  model  is  denoted  the 
non- equilibrium,  model.  A  second  and  a  relatively  simple  model  for  the  dust  phase  that  is 
incorporated  in  [19]  is  based  on  an  a-priori  assumption  of  velocity  equilibrium  among  the  two 
phases.  The  two-phase  mixture  is  described  by  the  incompressible  equations  and  an  additional 
equation  models  the  advection  of  the  dust  density  in  the  mixture  velocity  field,  i.e.,  a  tracer-like 
mechanism.  This  simple  model  is  denoted  the  equilibrium  model.  The  main  conclusion  of  the 
work  in  [19]  is  that  solutions  obtained  using  the  simplified  equilibrium  model  are  inherently 
different  from  solutions  obtained  using  the  non-equilibrium  model  described  above:  the  solutions 
exhibit  qualitatively  different  evolution  patterns  and  quantitatively  different  dust  density  fields. 
This  result  lead  us  to  focus  on  the  further  development  of  solvers  to  the  non-equilibrium  model 
only.  This  model  is  discussed  in  the  following  section. 
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2.2  Dusty  Gas  Model. 


For  the  dust  phase,  the  non-equilibrium  dusty  gas  model  of  Miura  and  Glass  [11]  is  used.  This 
model  assumes  that  the  dust  loading  is  sufficiently  small  such  that  the  volume  of  the  dust,  or 
suspension  phase,  can  be  neglected,  but  the  number  of  particles  is  large  enough  so  that  the 
suspension  phase  can  be  treated  as  a  continuum.  It  is  also  assumed  that  the  dust  does  not 
contribute  to  the  pressure  and  is  coupled  to  the  gas  phase  only  through  the  drag  terms.  Gravity 
forces  are  added  to  the  momentum  equations  since  these  are  important  for  many  applications. 
The  energy  conservation  equation  for  the  dust  phase  is  neglected  since  it  fully  decouples  from  the 
system  in  the  absence  of  a  gas  energy  equation  (and  hence  heat  transfer  terms)  and  it  does  not 
contribute  any  desirable  information  for  this  model.  The  incompressible  Euler  equations  for 
variable  density,  gas  phase  are  as  in  [6]  with  the  addition  of  the  appropriate  drag  terms  to  the 
momentum  equations.  The  system  of  equations  for  both  phases  based  on  these  assumptions  is: 


where, 


ut  +  (u  •  V)u 

Pt  +  (U  •  V)p 
Vu 
Ot  +  V  •  (crv) 
(erv)t  +  V  •  (dWT) 


(1) 


p  =  gas  density, 
u  =  {uj}  =  gas  velocity  vector, 
p  =  gas  pressure, 
a  =  dust  density, 
v  =  {^i}  =  dust  velocity  vector, 
g  =  gravity  vector, 

D  =  Di  =  single  particle  drag  vector, 
m  =  mass  of  a  single  dust  particle. 

The  particles  are  assumed  to  be  spheres  of  diameter  d.  In  past  work  we  used  the  drag  model  as 
suggested  in  Miura  and  Glass  [22].  That  is, 

D  =  £^2p(u  ~  v)  |  u  -  v  | CD,  (2) 

where  Cp  is  the  drag  coefficient  given  by 

CD  =  0.48  +  28i?e-0'85.  (3) 

Re  =  p  |  u  —  v  |  d/ p  is  the  Reynolds  number  and  p  is  the  gas  viscosity. 
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2.3  Extensions  to  Physical  Modeling. 


The  drag  model  given  by  Equations  (2-3)  is  based  on  Stokes  law  and  is  not  suitable  for  flows 
where  the  particles  experience  forces  that  accelerate  them  to  ultra-Stokesian  regimes.  For 
example,  setting  dust  particle  diameter  to  100  microns,  the  maximum  Reynolds  number  in  our 
sample  flows  ranges  in  order  of  magnitude  from  Remai  =  O(102)  (Gest  experiments)  to 
Remai  =  O(103)  (Priscilla  calculations),  which  are  2-3  orders  of  magnitude  (respectively)  larger 
than  the  values  for  which  Stokes  law  applies  (Re=  0(1)). 

We  have  implemented,  tested,  and  compared  results  for  three  additional  drag  models  in  this  work. 
Equation  (2)  above  is  denoted  Model  1  and  the  additional  models  chosen  are: 

(i)  the  model  suggested  by  Godoy  et  al.  [13]  (denoted  Model  2) 


where 


CD  = 


24 

Re 

10a 

Re 


Re  <  0.4207 
Re  >  0.4207 


s  = 


.133(log10  Re)  +  1.43 
.116(log10  Re)2  +  0.0544  log10 


0.4207  <  Re  <  1.9368 
Re  +  1.443Re  >  1.9368 


(ii)  the  model 


suggested  by  Wilson  [26]  (denoted  Model  3)  uses  the  drag  coefficient  given  by 


24 

Re 


24.  ,  4 

Re  fiFF 


Re  <  0.5 
Re  >  0.5 


(iii)  the  model  suggested  by  Chen  and  Pereira  [9]  (denoted  Model  4)  is  given  by, 

D=“^(u-v)fp 


(4) 


where  ap  is  the  particle  volume  concentration,  pp  is  the  mass  density  of  the  particle,  i.e.,  a  =  appp, 
fp  is  a  correction  factor  (equivalent  to  coefficient  of  drag),  depending  on  the  particle  Reynolds 
number;  r*  is  the  characteristic  response  time  of  the  particle  to  the  changes  in  fluid  motion,  given 

by  t*  =  and  fp  is  given  by  Boothroyd  (see  [9])  as  follows: 


1  +  0.15Re0687 
0.914Re°  282  +  0.0135Re 
0.0167Re 


0  <  Re  <  200 
200  <  Re  <  2500 
Re  >  2500 
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Section  3 
Numerical  Method 


3.1  General. 

The  solution  of  system  (1)  consists  of  coupling  two  distinctly  different  algorithms;  the  first  three 
equations  in  system  (1)  describe  the  evolution  of  a  variable  density,  gaseous  phase  of  the  flow  field 
under  the  constraint  of  divergence  free  velocity  field  (i.e.,  incompressible  flow).  The  solution  of 
these  equations  follow  closely  the  methodology  described  in  [6],  namely  the  projection  method. 
The  last  two  equations  in  (1)  describe  the  conservation  laws  for  the  dusty,  compressible  phase  of 
the  flow  field.  These  equations  are  solved  using  the  ideas  in  [12].  The  two  algorithms  are  coupled 
to  allow  the  appropriate  modeling  of  the  interaction  (or  coupling)  terms  of  the  equations,  namely 
the  drag  forces  D  =  D(u,  v,p).  Some  modifications  of  the  cited  algorithms  were  necessary  in 
order  to  construct  a  consistent  and  complete  methodology  for  the  non-equilibrium  model.  In  what 
follows  we  will  discuss  the  main  features  of  the  underlying  algorithms,  starting  with  a  brief 
description  of  the  code  of  [19]  and  following  with  the  necessary  extensions  applied  in  the  current 
work. 


3.2  The  Modified  Method  of  Krispin  and  Collins. 


The  basic  methodology  consists  of  first  solving  for  an  intermediate  velocity  field  u*  and  new 
density  field  using 


p^i/2g  _  vp»-i/ ^  _  lg"D"  +  g"+1D"+1' 

771  2 

(5) 

where  pn+L'z  in  Equation  (5)  is  the  average  of  pn  and  pn+1.  A  convection  algorithm  based  upon  a 
second-order  unsplit  Godunov  scheme  is  applied  to  extrapolate  velocities  to  cell  faces  at  tn+i/2.  At 
this  point  in  the  algorithm,  the  normal  velocities  on  cell  faces  are  centered  in  time  and  second 
order  accurate,  but  do  not,  in  general  satisfy  the  divergence  constraint.  In  order  to  enforce  the 
constraint  before  constructing  the  conservative  updates,  a  MAC  projection  [3]  is  applied  to  the 
velocities  at  the  intermediate  time  step  un+1/2  in  order  to  obtain  the  face-based  advection 
velocities  uADV .  Continuing  with  the  Godunov  scheme,  these  advection  velocities  are  then  used  to 
compute  the  conservative  updates,  namely  pn+1/2  and  the  velocity  field  u*  from  Eq.  (5)  which 
does  not,  in  general,  satisfy  the  divergence  constraint. 


ir  -  u" 
At 

pn+1 —  pn 


A  t 


+  [(u-  V)u 
+  [(u  •  V)p] 


ln+l/2  _ 


yi+ 1/2 


n+1/2  _ 


=  0, 


1  m  I  i  /o  . 


Then,  a  discrete  projection  P  is  applied  to  decompose  the  result  of  the  first  step  into  a  gradient  of 
a  scalar  potential  and  a  divergence-free  vector  field.  In  particular, 


un+1  —  u 
A t 

Vpn+1/2 

pn+l/2 


(6) 
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Given  some  initial  data  at  time  level  n  the  above  algorithm  requires  an  appropriate  approximation 
for  the  velocity  field  of  both  phases,  u  and  v,  and  the  densities  of  the  gas  and  dust,  p  and  <7, 
respectively,  at  time  level  n+1  for  the  calculation  of  the  source  term  crn+1Dn+1.  An  iterative 
procedure  was  constructed  that,  for  each  time  step,  starts  with  the  initial  guess  for  the  v  =  0 
iterate,  {a,p,  u,  v}n+1’"=0  =  {cr,  p,  u,v}n  and  repeats  the  computation  of  Equations  (3)  and  (4) 
either  until  some  convergence  criteria  is  met  or  for  a  sufficient,  fixed  number  of  iterations  u.  The 
update  for  the  gas  density,  pn+1,  is  given  at  the  end  of  the  calculation  in  Equation  (3)  whereas  the 
update  for  the  gas  velocity,  un+1,  is  the  end  result  of  the  projection  step  in  Equation  (4).  The 
update  for  the  velocity  field  of  the  dust  phase,  vn+1,  together  with  the  updated  dust  density,  an+1, 
are  computed  in  an  intermediate  step  which  is  embedded  in  the  above  iterative  procedure.  This 
step  is  based  on  the  ideas  in  [12],  as  mentioned  above,  and  it  is  discussed  in  [19]  in  great  detail. 
The  main  limitation  of  the  dust  phase  integration  procedure  is  the  fact  that  it  is  using  the 
operator-split  technique  both  for  spatial  and  source  integrations.  This  technique  was  modified  as 
discussed  in  the  following  section. 


3.3  Algorithmic  Extensions. 

We  started  this  work  by  incorporating  Colella’s  unsplit  methodology  [11]  in  the  dust  phase  of  our 
original  code  [19].  As  a  consequence,  source  integration  in  the  dust  phase  can  be  performed  in  an 
unsplit  fashion  which  is  a  necessary  requirement  for  stable  and  accurate  integration  of  stiff  source 
terms  (see,  e.g.,  conclusions  of  the  studies  in  [18]  and  [20]). 

The  new  conservative-difference  corrector  step  of  the  dust  phase,  written  in  two-space  dimensions, 
is  given  by, 


WM  - 


1/2  -[^71+1/2 

.j+1/2,*  r  j-l/2,fc 


CX.J, 


r^n-f  1/2 

IV  j, *+1/2 


-  G j,k-  l/2"+1/2]  +  LS(W”+>,  WJt) 


(7) 


where  W,  F,  G,  and  S  are  defined  in  the  obvious  way  (using  Equation(l)),  ax  =  ay  —  At, 
and  L  denotes  a  source-term  integration  operator.  We  tested  two  operators,  given  by 


£.  =  Y(S(ws>)  +  s(wy) 


(8) 


and 

L2  =  AtS(WJ+1). 


(9) 


Other  algorithmic  extensions  involve  cycling  and  subcycling  needed  to  iterate  on  the  implicit  (or 
semi-implicit)  variables  that  are  coupled  between  the  two  phases,  in  particular  the  velocity  fields. 
Several  user-defined  iteration  counts  are  built  into  the  code  and  they  are  problem  dependent. 
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Section  4 
Numerical  Results 


4.1  Convergence  Study. 

The  rate  of  convergence  for  the  method  discussed  above  is  investigated  for  a  smooth  initial  stream 
function  and  density  field  inside  a  unit  square  defined  exactly  as  in  [6], 

0,  y )  =  n~1sin2(iry)sin2(Trx),  p(y)  =  1  -  \tanh(y  -  -).  (10) 

4  2 

The  initial  velocities  are  defined  in  the  usual  way  by  u  =  dijj/dy  and  v  =  -dijj/dx.  The  initial 
dust  density  a  =  O.Olp,  and  the  initial  dust  velocity  v  =  u.  The  velocities  for  both  phases  satisfy 
homogeneous  Dirichlet  boundary  conditions  and  the  density  for  both  phases  is  determined 
through  a  second-order  extrapolation  from  the  outermost  interior  cells.  The  gravity  term  was 
removed  for  the  purposes  of  this  study.  The  solution  is  computed  on  uniform  grids 
Aa;  =  A y  =  l/2n  for  2n  uniform  time  steps  at  a  CFL  =  0.5  for  n  =  4,  5,  6,  7.  The  second  normal 
residual  is  measured  between  grids  of  adjacent  resolution  and  is  used  to  calculate  the  rate  of 
convergence  which  are  summarized  in  Table  1.  The  gas  density  and  velocity  fields  for  the  coupled 
system  retains  a  second-ordered  accuracy  with  a  slight  deterioration  similar  to  that  reported  in  [6] 
for  the  single-phase  system.  The  convergence  rate  for  the  particle  density  and  velocity  fields 
indicate  an  accuracy  between  first  and  second  order.  This  is  due  to  the  difficulties  of  accurately- 
predicting  the  flux  terms  in  the  particle-phase,  particularly  in  regions  where  vacuum  states  or 
colliding  particle  streams  occur.  Such  limitations  of  the  particle-phase  model  are  also  discussed  in 
[12]  for  compressible  flows.  Modifications  to  reduce  the  limitations  of  the  particle-phase  model  are 
currently  under  investigation 


Table  1.  Rates  of  convergence. 


Variable 

32-64 

Rate 

64-128 

Rate 

128-256 

P 

5.612e-04 

2.12 

1.245e-04 

1.94 

2.217e-05 

u 

3.025e-03 

1.87 

8.628e-04 

1.77 

2.746e-04 

a 

5.908e-04 

1.34 

3.271e-04 

1.34 

1.810e-04 

V 

5.589e-03 

1.42 

2.791e-03 

1.37 

1.477e-03 

4.2  Validation. 

In  the  same  manner  as  that  of  Krispin  and  Collins  [19],  the  numerical  scheme  with  the  algorithmic 
extensions  described  in  section  3.3  is  initially  tested  by  simulating  low  density  bubbles  containing 
a  dust  phase.  The  gas  density  of  an  initially  spherical  bubble  is  specified  to  be  one  tenth  the 
density  of  the  surrounding  ambient  gas.  The  density  of  the  dust  particle  phase  is  specified  to  be 
one  hundredth  of  the  density  of  the  bubble’s  gas  phase.  The  bubble  is  initially  at  rest. 
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r  <m) 


(a)  0  seconds 


(b)  2.5  seconds 


(c)  5.0  seconds 


Figure  1.  Gas  and  dust  density  contours  for  simulation  of  a  dusty  bubble  with  100  micron  dust 
particles. 


Figure  1  shows  the  result  for  an  axisymmetric  calculation  for  a  dusty  bubble  with  a  100  micron 
particle-phase.  The  simulation  is  from  time  zero  to  5.0  seconds.  The  color  scale  depicts  the  density 
of  the  gas-phase,  where  blue  indicates  lower  densities  and  red  higher  densities.  Ten  contour  lines 
map  out  the  loglO  density  of  the  dust-phase.  The  velocity  vectors  depict  the  gas-phase  velocities. 

Gravity  is  oriented  in  the  usual  direction.  As  the  low  density  gas  bubble  rises,  the  dust  near  the 
top  of  the  bubble  rises  as  well,  entrained  in  the  resulting  gas  velocity  field.  The  dust  in  the  region 
around  the  lower  portion  of  the  bubble  falls  due  to  the  gravity  field.  As  expected  rate  of  fallout  at 
of  the  dust  in  the  lower  portion  of  the  bubble  is  found  to  be  sensitive  to  the  specified  particle 
diameter.  The  test  case  in  Figure  1  has  a  grid  resolution  of  32x64  for  a  physical  domain  of  10x20 
meters. 

Figure  2  shows  a  three-dimensional  dusty  bubble  simulation  where  in  this  case,  the  particle-phase 
is  modeled  as  an  advected  scalar  in  the  same  manner  as  the  AD  scheme  described  in  Krispin  and 
Collins  [19]  (although  here  we  now  utilize  the  projection  method  described  in  section  3.3).  The 
lefthand  frame  depicts  a  blue  gas  density  isosurfaces  and  black  log  particle  density  isosurfaces. 

The  righthand  frame  shows  a  slice  of  the  flow  field  in  the  x-z  plane  at  approximately  y  =  10 
meters.  Similar  to  the  two-dimensional  simulations,  the  color  scale  depicts  the  gas  density,  the 
black  line  contours  map  the  log  of  the  dust  density,  and  the  velocity  vectors  indicate  the  motion  of 
the  gas  phase.  The  gas  density  in  the  bubble  is  less  than  that  of  the  surrounding  ambient  gas. 

This  causes  the  bubble  to  rise  while  continuously  entraining  the  heavier  ambient  gas. 

In  the  same  manner  as  the  two-dimensional  axisymmetric  cases,  the  simulation  is  from  time  zero 
to  5.0  seconds.  The  dust  density  for  infinitely  small  particles  is  computed  through  a  single 
advection  equation,  i.e.,  the  dust  particles  are  completely  entrained  into  the  gas-phase  flow  and 
are  simply  advected  with  the  gas. 

This  test  case  has  a  coarse  grid  of  32x32x32  cells  for  a  domain  of  20  meters  cubed.  A  runtime  of 
approximately  0.5  CPU  hours  on  a  400  MHz  Pentium  II  processor  is  required  to  advance  the 
solution  to  5.0  seconds. 
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In  summary,  both  the  multi-phase  model  and  single  phase  advection  model  result  in  the  same 
qualitative  trends  for  rising  dusty  bubbles  as  those  described  in  Krispin  and  Collins  [19]  for  both 
axisymmetric  and  three-dimensional  simulations. 


4.3  Confined  Explosions. 

For  a  number  of  reasons  including  the  availability  of  model  input  data  and  experimental  results, 
the  bulk  of  our  effort  was  spent  on  the  simulation  of  confined  explosions.  The  problem  however, 
involves  both  compressible  and  incompressible  flow  regimes,  so  the  simulation  was  undertaken  in 
sections.  Initially,  the  compressible  portion  of  the  simulation  was  accomplished  by  simply  using 
the  plume  size  and  mass  distribution  prediction  by  Applied  Research  Associates’  (ARA) 
STructural  Expulsion  and  Plume  (STEP  2.3)  model.  [14]  Using  STEP  data  as  a  starting  condition 
for  the  incompressible  model  proved  to  be  a  poor  choice  for  for  initial  conditions  however.  The 
plume  predicted  by  STEP  as  early  as  8.4  seconds  was  negatively  buoyant  and  fell  quickly  to  earth 
after  mapping  to  CATAPULT.  Figure  3  shows  the  evolution  of  the  plume  when  mapped  from 
STEP  input  from  8.4  seconds.  One  of  the  primary  reasons  for  this  collapse  is  most  likely  due  to 
the  lack  of  a  stratified  atmosphere  in  the  STEP  code  which  models  the  atmosphere  as  an  infinite 
room  with  a  constant  air  density.  STEP  is  also  based  on  empirical  data  and  linearized  equations 
which  result  in  inaccurate  predictions  of  dust  and  gas  distribution  inside  the  plume. 

After  determining  that  STEP  provided  a  poor  density  distribution  for  the  input  to  the 
incompressible  model,  it  was  decided  that  a  realistic  simulation  of  the  plume  required  modeling 
the  compressible  phase  of  the  problem  as  well  as  the  later  incompressible  stages.  The  compressible 
phase  of  the  plume  was  modeled  using  vent  data  that  was  also  provided  by  STEP.  However,  since 
it  is  more  difficult  to  predict  realistic  3-D  velocity  and  density  distributions  than  it  is  to  predict 
velocity  and  density  at  a  point  (the  vent),  it  was  assumed  that  the  vent  data  would  be  more 
accurate  than  the  distribution  predictions  used  earlier. 

The  results  from  the  DIPOLE  EAST  experiment  in  New  Mexico  indicated  that  the  plume  should 
rise  to  approximately  120  meters  in  height  after  10  seconds,  while  spreading  to  a  width  of 
approximately  35  meters  in  the  same  time.  The  results  then  show  that  the  plume  rises  to  a  height 
of  approximately  325  meters  and  an  average  width  of  180  meters  over  the  following  couple  of 
minutes  before  dissipating.  STEP  predictions,  which  are  based  primarily  upon  empirical  data  fits 
from  similar  experiments  also  predict  that  the  plume  initially  rises  quickly  (0-10  seconds)  to  a 
level  of  approximately  150  meters  and  spreads  to  a  diameter  of  20  meters  over  the  same  time. 
These  results  correspond  to  spreading  rates  ( drjdz )  for  the  plume  of  0.146  and  0.073  for  the 
experiment  and  STEP  respectively  over  the  first  ten  seconds  of  the  flow. 


4.3.1  Compressible  Phase. 

Initially,  an  axisymmetric  compressible  model  was  used  to  simulate  the  compressible  portion  of 
the  experiment,  but  problems  were  soon  discovered  with  that  approach.  Using  both  Logicon  R&D 
Associates  report  on  the  DIPOLE  EAST  159  experiment[2]  and  plume  height  predictions  from 
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Figure  3.  Gas  density  color  contours  with  Log(dust  density)  contours  overlain  for  incompressible 
run  initialized  with  STEP  data. 


11 


STEP  as  a  basis,  it  quickly  became  apparent  that  the  axisymmetric  model  was  predicting  plume 
heights  that  were  entirely  too  high. 

Several  axisymmetric  runs  using  full  AMR  were  undertaken  and  one  was  conducted  out  to  more 
than  two  seconds.  In  every  axisymmetric  case,  the  plume  spreading  rate  was  extremely  small  and 
the  plume  height  predictions  were  advancing  too  rapidly.  Perturbations  were  introduced  to  both 
the  ambient  atmosphere  and  radial  velocity  but  did  not  significantly  affect  this  outcome.  Despite 
the  efforts  to  induce  instability  in  the  flow  and  enhance  plume  spreading  and  slowing,  all 
axisymmetric  compressible  runs  resulted  in  unrealistically  high  and  narrow  plumes. 

Figure  4  shows  the  plume  heights  predicted  by  axisymmetric,  and  three  dimensional  simulations 
as  well  as  the  baseline  prediction  by  STEP.  The  curves  clearly  show  that  axisymmetric 
simulations  even  when  performed  at  high  resolution  with  AMR  are  not  able  to  capture  the  physics 
of  the  problem,  whereas  the  three  dimensional  simulation  is  better  able  to  do  so  even  when  run  on 
relatively  coarse  resolution  grids.  The  axisymmetric  runs  displayed  a  large  variety  of  vortex  rings 
and  roll-ups  along  an  otherwise  very  narrow  jet-like  structure.  The  vortex  formations  along  the  jet 
flow  can  be  seen  in  Figure  5  which  shows  gas  density  for  the  entire  plume  at  several  times  in  one 
simulation  as  well  as  a  closer  view  of  the  base  of  the  plume  at  0.326  seconds.  Clearly  evident  in 
the  figure,  the  plume  exhibits  little  spreading  and  the  plume  height  is  considerably  higher  than 
that  seen  in  either  experimental  results  or  STEP  predictions  shown  in  Figure  4.  The  fourth  panel 
in  Figure  5  exhibits  the  considerable  turbulence  that  is  evident  along  the  centerline  of  the  domain 
and  is  found  at  all  heights  along  the  jet.  Despite  this  turbulent  energy  sink,  the  plume  did  not 
slow  or  spread  appreciably,  indicating  that  the  main  mode  for  energy  transfer  occurs  elsewhere  in 
the  problem. 

Analysis  of  the  plume  spreading  rate  in  the  early  times  of  the  axisymmetric  simulations  shows 
that  it  is  between  0.1  and  0.3  (depending  on  where  the  edge  of  the  plume  is  determined  to  be) 
which  is  right  in  line  with  the  spreading  rates  in  the  near  field  reported  by  Kouros  et  al.[16]  and 
the  experimental  results.  The  plume  tip  is  also  only  at  35  meters  or  58  d0  downstream  at  this 
point.  The  problem  arises  at  later  times  when  the  plume  has  not  continued  to  spread  away  from 
the  axis.  At  approximately  a  second  after  the  blast,  the  plume  has  continued  to  rise  with  little  or 
no  spreading,  reducing  the  spreading  rate  to  0.032.  The  tip  of  the  plume  at  this  point  is  at  nearly 
200  meters  or  316  d0  downstream.  Experimental  jets  and  plumes[21,  16]  generally  show  that 
axisymmetry  begins  to  break  between  50  and  100  d0  and  a  three  dimensional  simulation  is 
required  to  capture  the  correct  asymmetric  motions. 

Since  attempts  to  induce  instability  in  the  axisymmetric  simulations  by  introducing  density  and 
inflow  velocity  perturbations  were  unsuccessful,  it  became  apparent  that  the  axisymmetric 
formulation  itself  was  not  capturing  the  full  physics  in  the  problem  and  full  three  dimensional 
simulations  were  undertaken.  The  3D  simulation  was  run  to  a  time  of  10  seconds  before  being 
mapped  to  an  incompressible  grid.  A  series  of  plots  show  gas  and  dust  density  isosurfaces  for  10 
times  throughout  the  3-D  simulation  in  Figure  6.  The  axisymmetry  in  the  problem  is  already 
beginning  to  break  in  the  first  frame  at  a  time  of  0.73  seconds  and  the  flow  is  definitely 
asymmetric  by  frame  2  which  is  at  a  time  of  1.4  seconds.  This  trend  toward  asymmetry  is 
observed  in  both  the  dust  and  gas  isosurfaces  throughout  the  simulation  as  the  plume  height 
reaches  more  than  120  meters  by  a  time  of  8.1  seconds,  matching  STEP  and  DIPOLE  East  results 
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Time  (seconds) 


Figure  4.  Plume  heights  from  axisymmetric,  3D,  and  STEP  simulations  (above)  and  a  blow  up 
view  of  axisymmetric  and  3D  simulations  at  early  times  (below). 
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(a)  Plume  density  at  time  0.190  sec¬ 
onds. 
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(b)  Plume  density  at  time  0.326  sec¬ 
onds. 
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(c)  Plume  density  at  time  1.021  sec¬ 
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Figure  5.  Plume  gas  density  at  various  times  from  axisymmetric  simulations 
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(a)  Time  =  0.73  seconds 


(b)  Time  =  1.4  seconds 


(c)  Time  =  2.1  seconds 


(d)  Time  =  2.8  seconds 


(e)  Time  =  3.6  seconds 


(f)  Time  =  4.4  seconds 


(g)  Time  =  5.1  seconds 


(h)  Time  =  5.9  seconds 


(i)  Time  =  6.9  seconds  0)  Time  =  8.1  seconds 

Figure  6.  Gas  density  isosurface  (left  panels)  colored  according  to  the  Log  of  dust  density  and  dust 
density  isosurface  (right  panels)  colored  according  to  gas  density. 
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quite  well.  Three  dimensional  results  began  diverging  noticeably  from  axisymmetry  at  around  half 
a  second  after  the  blast,  and  this  divergence  significantly  affected  the  plume  height  and  spreading 
rate.  As  is  evident  in  Figure  4,  both  axisymmetric  and  3D  simulations  agree  very  closely  for  the 
first  0.33  seconds  of  the  flow.  The  spreading  rate  for  both  simulations  are  likewise  similar  up  to 
this  time.  At  around  0.33  seconds  however,  the  tip  of  plume  in  the  3D  simulation  begins  to 
expand  and  move  off  the  centerline  of  the  domain  which  appears  to  be  the  primary  cause  of  the 
slowing  in  the  3D  simulation  and  in  real  plumes.  It  should  be  noted  that  the  divergence  does  not 
occur  until  the  tip  of  the  plume  has  achieved  a  height  of  nearly  lOOdo,  which  is  as  much  as  most 
jet  flow  simulations  and  experiments  in  the  literature  capture  of  the  flow  at  all. 


Non-dimensionsal  Time 


Figure  7.  Plume  heights  vs.  square  root  of  non-dimensional  time  ( Vjtd~ 1)0-5. 

Johari  et  al.  [15]  have  shown  that  impulsively  started  jets  display  a  starting  vortex  that  scales  in  a 
self  similar  manner  as  the  square  root  of  non-dimensional  time  (fVjd^1).  A  plot  of  the  data  shown 
in  Figure  7  from  both  axisymmetric  and  fully  3D  simulations  shows  that  all  initially  display  a 
linear  trend  similar  to  that  seen  by  Johari  et  al.  The  slope  of  the  axisymmetric  simulations 
diverges  quickly  from  the  STEP  prediction  and  the  3D  simulation,  but  all  three  display  a  mostly 
linear  trend  through  the  majority  of  the  simulation.  As  the  vent  inflow  velocity  drops  rapidly  to 
zero  a  reversal  of  non-dimensional  time  occurs  and  is  responsible  for  the  switch-backs  seen  in  the 
STEP  and  3D  curves.  Both  the  STEP  and  3D  curves  fall  slightly  below  the  results  reported  by 
Johari  et  al.  (also  shown  on  the  figure),  while  the  axisymmetric  run  falls  above  all  the  others  by  a 
considerable  margin.  The  slope  reported  by  Johari  et  al.  is  approximately  2.7  whereas  the  3D 
simulation  predicts  a  slope  of  roughly  1.95.  The  differences  in  slopes  may  indicate  that  plume 
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penetration  slows  at  higher  Reynolds  numbers  or  that  compressible  effects  and  gravity  are  slowing 
the  plume. 

Typical  averaged  velocity  profiles  for  a  steady  state  turbulent  jet  follow  a  self-similar  solution  that 
has  been  fairly  well  established.  [25]  When  scaled  by  the  axial  velocity  (W/Wax)  and  plotted 
against  axial  distance  over  radial  distance  (Z/R),  the  velocity  profile  collapses  to  sech2(Z/R). 
Atassi  et  al.  [4]  have  shown  that  transient  jets  can  be  expected  to  follow  a  modified  self-similarity 
as  well.  It  has  been  shown  that  when  far  enough  removed  from  the  jet  source,  the  velocity  profile 
of  an  unsteady  jet  is  the  same  as  that  seen  in  a  steady  jet.  [16]  The  assumption  however,  is  that 
the  jet  oscillates  around  some  non-zero  average,  which  is  somewhat  different  than  the  case 
presented  here. 

Since  the  velocity  profile  is  scaled  by  the  velocity  along  the  axis,  it  is  not  unreasonable  to  expect 
an  unsteady  decaying  plume  to  behave  similarly.  Figure  8  shows  the  self-similar  velocity  structure 
from  the  3D  compressible  simulation  of  the  plume  over  several  different  time  steps.  There  is  a 
significant  amount  of  scatter  in  the  data,  but  the  velocity  generally  follows  the  expected  velocity 
profile  which  is  plotted  as  a  solid  line.  Initially,  the  profiles  from  higher  altitudes,  which  are 
plotted  as  stars,  are  scattered  fairly  evenly  above  and  below  the  theoretical  profile,  but  at  later 
times,  they  fall  consistently  below  the  theoretical  curve,  especially  at  larger  values  of  R/Z.  This  is 
most  likely  indicative  of  gravity  beginning  to  effect  the  flow.  The  tightening  of  the  profile  at  higher 
altitudes  also  lends  credence  to  the  idea  that  instability  increases  away  from  the  axis  of  symmetry. 
The  general  conformance  to  the  self-similar  profiles  for  steady  jets  indicates  that  the  overall  flow 
still  conforms  to  an  ensemble  axial  symmetry.  However,  the  comparisons  of  3D  and  axisymmetric 
simulations  with  empirical  models  and  experimental  results  show  undeniably  that  the  strict 
enforcement  of  axial  symmetry  produces  incorrect  values  for  both  plume  height  and  spreading 
rate.  This  leads  us  to  the  conclusion  that  compressible  axisymmetric  simulations  are  inadequate  in 
capturing  the  real  physics  in  jet  and  plume  problems.  While  results  may  be  fairly  accurate  in  the 
near  field,  by  the  time  the  flow  has  reached  50  do,  the  axial  symmetry  is  no  longer  appropriate. 

Volume  entrainment  in  a  plume  or  jet  can  be  computed  as  the  volume  of  gas  that  is  within  the 
plume  flow  field.  This  reduces  to  a  summation  of  gas  density  multiplied  by  vertical  velocity  and 
grid  cross  sectional  area  given  as 

nx  ny 

Q  =  £  5Z  Pijwijdxdy  (11) 

i=l j= 1 

and  is  computed  at  every  z  level  in  the  flow.  This  yields  a  volume  flow  rate  at  each  level  in  the  z 
direction  at  any  given  time.  Bremhorst  and  Hollis  [7]  have  shown  that  the  entrainment  rate  for  a 
fully  pulsed  jet  is  nearly  twice  that  observed  in  steady  jets  in  the  near  field.  They  also  found  that 
the  volume  entrainment  rate  scaled  as  a  nearly  linear  function  of  x/d.  Figure  9  is  a  contour  plot  of 
volume  entrainment  as  a  function  of  time  and  distance.  The  volume  entrainment  during  the 
compressible  phase  of  the  3D  plume  simulation  is  roughly  linear  at  lower  values  of  z/d  as  noted  by 
the  regular  spacing  in  contour  levels  up  to  a  level  of  z/d  «  75.  The  entrainment  curve  exhibits  a 
strongly  increasing  trend  at  higher  levels  of  z/d  before  dropping  precipitously  at  the  top  of  the 
plume.  The  entrainment  rate  is  also  a  function  of  time  as  well.  As  time  progresses,  the 
entrainment  increases  in  volume  as  the  main  area  of  entrainment  propagates  away  from  the  vent. 
The  near  field  volume  entrainment  rate  is  higher  than  that  reported  for  steady  jets  at  early  times, 
but  less  than  the  rate  reported  by  Bremhorst  and  Hollis,  falling  at  approximately  1.5  times  the 


17 


0.00  0.05  0.10  0.15  OJO 

R/Z 
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(c)  Self-similar  velocity  profile  at 
5.93  seconds. 
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(b)  Self-similar  velocity  profile  at 
4.35  seconds. 
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(d)  Self-similar  velocity  profile  at 
8.18  seconds. 


Figure  8.  Self-similar  plume  velocity  profiles  from  various  times.  Stars  indicate  velocities  from 
higher  altitudes  while  diamonds  represent  lower  altitudes  in  the  plume.  The  line  is  the 
theoretical  jet  velocity  profile. 
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steady  rate.  As  the  entrainment  area  propagates  away  from  the  vent  however,  the  rate  of 
entrainment  increases  non-linearly  to  a  maximum  value  of  114  Q/QE  at  154.7  Z/dQ,  where  Q/QE 
is  the  volume  flow  rate  normalized  by  the  inflow  rate.  Although  Bremhorst  and  Hollis  did  not 
report  volume  rates  at  these  distances,  if  their  reported  entrainment  rate  continues  in  a  linear 
fashion,  it  would  be  approximately  117  Q/QE  at  154.7  x/d0.  The  increasing  entrainment  rate  in 
the  far  field  may  suggest  that  the  optimal  mixing  region  in  conventional  jets  occurs  at  distances 
beyond  those  investigated  in  previous  experiments. 

Longmire  and  Eaton[21]  have  shown  that  for  two-phase  flows,  particles  congregate  in  areas  of  gas 
divergence  while  they  are  likely  to  be  sparse  in  areas  of  convergence.  The  largest  areas  of 
divergence  in  this  flow  also  have  a  tendency  to  show  higher  concentrations  of  dust.  Figure  10 
shows  an  Y-Z  cross  section  of  the  plume  at  times  of  5.1  and  6.9  seconds  with  the  shaded  contours 
showing  gas  divergence  and  the  line  contours,  dust  density.  The  correlation  is  not  exact,  but 
higher  dust  densities  are  generally  in  the  darker  or  divergent  regions  of  the  plot  just  as  predicted 
by  Longmire  and  Eaton. 
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Figure  10.  Gas  divergence  in  shaded  contours  with  dust  contour  lines  overlaid. 

4.3.2  Incompressible  Phase. 

Once  the  compressible  portion  of  the  flow  was  concluded,  the  3D  results  were  mapped  to  both 
axisymmetric  and  3D  incompressible  domains.  Although  axisymmetry  has  proven  to  be 
inadequate  in  the  compressible  region,  it  was  not  clear  that  this  would  also  be  the  case  for  the 
later  stages  of  the  flow.  Since  the  inflow  to  the  domain  from  the  vent  ended  at  roughly  the  6 
second  mark  of  the  compressible  phase  and  gas  velocities  have  dropped  considerably,  the  problem 
then  becomes  an  incompressible,  buoyant,  stratified,  multi-phase  flow. 

There  is  some  question  as  to  where  the  flow  should  be  mapped  from  compressible  to 
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incompressible  regimes.  For  steady  flows,  the  general  rule  for  transition  to  incompressible  flow  is 
that  M  «  1.  The  flow  drops  below  a  Mach  number  of  0.1  at  a  time  of  approximately  6  seconds, 
but  the  forcing  continues  until  shortly  after  6  seconds.  In  order  to  handle  the  inflow  in  a 
consistent  manner,  it  was  thought  best  to  model  the  problem  as  compressible  at  least  until  the 
forcing  from  the  vent  was  finished.  The  compressible  portion  of  the  simulation  was  run  out  to  a 
time  of  10  seconds  and  it  was  mapped  to  the  incompressible  regime  at  times  of  6.9  seconds,  8.1 
seconds  and  10  seconds  for  comparison.  Although  the  solution  was  not  quite  divergence  free  at  the 
time  of  the  mapping,  the  compressible  solution  was  sufficiently  convergent  that  the  input  did  not 
destabilize  the  incompressible  code. 

Axisymmetric  runs  were  made  by  averaging  the  3D  compressible  results  axially  and  mapping 
them  into  a  domain  that  was  75  by  450  meters  in  size  using  resolutions  of  both  16  by  128  and  32 
by  256.  Since  the  domain  is  axisymmetric,  the  effective  size  of  the  simulation  is  150  by  450 
meters.  The  3D  runs  utilized  the  same  effective  domain  size  as  the  axisymmetric  runs  and  a 
resolution  of  32  by  32  by  128  grid  points. 

All  simulations  were  computed  out  to  a  time  of  125  seconds  after  the  initial  blast,  and  the  results 
were  somewhat  surprising.  Plume  height  in  the  axisymmetric  simulations  was  a  strong  function  of 
mapping  time,  varying  from  nearly  400  meters  when  mapped  at  6.9  seconds  into  the  high 
resolution  axisymmetric  model  to  270  meters  when  mapped  from  a  time  of  10  seconds.  The 
axisymmetric  model  consistently  predicted  higher  plume  heights  than  did  the  fully  3D  simulations 
started  from  any  time.  At  higher  resolutions,  the  axisymmetric  results  were  even  further  from  the 
3D  predictions.  The  reason  for  this  appears  to  be  the  ability  to  resolve  a  tighter  vortex  ring  at  the 
top  of  the  plume  in  the  higher  resolution  cases.  Higher  resolution  runs  predicted  a  higher  thinner 
plume  than  did  the  coarser  runs  from  the  same  time,  which  seems  to  be  entirely  due  to  vortex  ring 
at  the  top  of  the  plume.  While  vortex  rings  are  present  in  all  the  axisymmetric  simulations,  they 
exist  only  in  a  loose  manner  in  the  3D  simulations  and  break  up  before  they  can  propagate  very 
far.  Figure  11  shows  the  predicted  plume  heights  from  high  and  low  resolution  axisymmetric  runs 
as  well  as  those  from  3D  simulations  started  from  the  three  different  initial  conditions.  Note  that 
the  high  and  low  resolution  axisymmetric  runs  match  very  well  at  early  times  before  the  vortex 
ring  breaks  up  in  the  coarse  resolution  simulations.  Also,  while  the  3D  simulations  do  not  vary 
greatly  in  predicted  plume  height,  but  the  axisymmetric  simulations  are  highly  varied.  The 
experimental  results  suggest  that  the  plume  should  rise  to  a  height  of  approximately  290  meters 
which  is  best  matched  by  the  axisymmetric  runs  mapped  at  a  time  of  8.1  seconds.  Since  the  other 
axisymmetric  results  are  so  erratic,  this  is  most  likely  just  a  coincidence.  The  large  envelope  of 
plume  height  predictions  by  the  axisymmetric  model  suggest  that  great  care  should  be  used  before 
settling  for  an  axisymmetric  simplification  to  buoyant  plume  flows.  The  relative  consistency  seen 
in  the  3D  simulations  also  suggests  that  the  exact  map  over  point  from  compressible  to 
incompressible  flow  is  not  highly  sensitive  if  a  full  3D  simulation  is  utilized.  At  least  some  of  the 
difference  between  the  height  prediction  by  the  3D  simulation  and  the  experimental  results  is 
likely  a  result  of  inflow  conditions  in  the  compressible  model.  Since  these  conditions  were  provided 
by  another  model  and  not  actual  measurements,  an  agreement  in  heights  between  experiment  and 
3D  simulations  of  less  than  13%  was  quite  acceptable. 

While  the  axisymmetric  assumption  proved  to  be  less  disastrous  to  the  solution  in  this  phase  of 
the  problem  than  it  was  in  the  compressible  phase,  it  is  still  not  a  good  assumption  in  this  flow.  A 
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side  by  side  comparison  of  frames  from  the  low  resolution  axisymmetric  run  and  the  3D  run 
mapped  from  6.9  seconds  shown  in  Figure  12  exhibits  the  difference  in  distribution  between  the 
two  simulations.  The  axisymmetric  constraint  forces  the  flow  to  move  in  the  form  of  a  vortices 
directly  atop  a  steady  plume,  whereas  the  unconstrained  3D  flow  shows  a  vortical  structure  atop  a 
generally  chaotic  velocity  field  within  the  plume.  A  series  of  plots  from  the  fully  3-D  simulation  is 
show  in  Figure  14  and  shows  a  very  realistic  evolution  of  the  dust  density  over  the  incompressible 
life  of  the  plume,  spanning  from  the  original  mapping  to  stagnation  and  dust  fall  out..  In 
addition,  the  velocity  distributions  in  the  two  simulations  are  quite  dissimilar  and  lead  to 
differences  in  dust  distribution  as  well.  The  axisymmetric  simulations  tend  to  loft  more  dust  to 
higher  altitudes  than  do  the  3D  simulations  keeping  the  dust  closer  to  the  centerline  of  the 
domain  than  is  seen  in  3D.  Figure  13  shows  the  vertical  and  horizontal  dust  distribution  in 
kilograms  for  the  3D  and  coarse  axisymmetric  models  started  from  a  time  of  6.9  seconds  and 
highlights  this  tendency.  In  particular,  note  the  spike  at  the  top  of  the  axisymmetric  plot  of 
vertical  distribution  that  shows  a  disproportional  amount  of  dust  in  the  tip  of  the  plume. 

Through  comparison  with  experimental  results  and  empirical  models,  it  has  been  shown  that 
plume  phenomenon  are  not  well  reproduced  by  axisymmetric  compressible  models.  Even  with 
extremely  high  resolution,  the  basic  physics  of  the  problem  is  not  represented  by  an  axisymmetric 
assumption.  Although  the  axisymmetric  and  3D  results  agreed  well  in  the  near  field  region,  in  the 
far  field  where  the  plume  tip  was  located  after  less  than  a  second,  only  the  3D  simulation  was  able 
to  accurately  reproduce  experimental  results.  The  plume  height  and  spreading  rate  were  both 
captured  in  a  much  better  manner  by  a  fully  three  dimensional  model  at  later  stages  of  the  flow  in 
the  compressible  regime.  The  ensemble  velocity  profiles  of  the  3D  simulation  fit  the  expected 
theoretical  profiles  to  a  reasonable  extent,  but  the  actual  asymmetry  in  the  flow  proved  to  be  a 
major  contributor  to  the  resulting  plume  height. 

The  entrainment  rates  observed  in  the  compressible  phase  of  the  3D  flow  were  slightly  lower  than 
the  rates  reported  by  Bremhorst  and  Hollis  [7]  for  near  field  values,  but  began  to  increase  quickly 
in  the  far  field.  The  increasing  entrainment  rate  in  the  far  field  might  suggest  that  increased 
mixing  could  be  achieved  in  conventional  jet  flows  by  spreading  the  same  volume  flow  through  a 
number  of  smaller  jets. 

Finally,  the  incompressible  phase  of  the  flow  has  also  proven  to  be  poorly  suited  to  the 
axisymmetric  constraint.  The  axi-symmetry  in  the  problem  continues  to  predict  large  scale 
coherent  structures  in  the  flow,  that  are  most  likely  not  present.  In  addition  to  incorrectly 
capturing  the  velocity  fields,  this  tendency  toward  coherent  structures  in  the  axisymmetric  model 
can  significantly  affect  the  distribution  of  dust  in  a  multi-phase  flow.  Plume  height  and  general 
shape  were  also  found  to  be  inadequately  captured  with  an  axisymmetric  model,  especially  when 
run  at  higher  resolutions. 
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(a)  Time  after  blast  is  40  seconds. 
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(b)  Time  after  blast  is  90  seconds. 

Figure  12.  Gas  density  in  shaded  contours  with  an  isosurface  of  dust  density  in  translucent  blue 
shown  from  two  times  during  the  simulation.  The  3D  results  are  shown  on  the  left 
and  the  axisymmetric  results  displayed  in  three  dimensional  form  are  on  the  right. 
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(a)  3D  dust  distributions 


(b)  Axisymmetric  dust  distributions 


Figure  13.  Horizontal  (above)  and  vertical  (below)  dust  distribution  in  3D  and  axisymmetric 
simulations  started  from  6.9  seconds  after  the  blast.  These  distributions  are  from 
a  time  of  30.8  seconds  after  the  blast. 

4.4  Unconfined  Explosions. 

4.4.1  Introduction. 

There  was  very  little  existing  data  to  use  for  comparison  to  simulated  results  presented  here. 
STEP  provides  an  unconfined  explosion  model  component,  but  it  is  relatively  undeveloped  and 
was  useful  only  to  provide  initial  conditions  to  CATAPULT.  The  problem  was  mapped  into  the 
incompressible  component  of  C ATAP ULT  from  the  STEP  results  using  a  simulated  dust  and  gas 
distribution.  STEP  output  contained  only  a  constant  density  component  for  the  gas  and  no  dust 
component,  so  gas  and  dust  distributions  similar  to  those  in  hemispherical  top  of  a  plume 
predicted  by  STEP  was  utilized  as  the  initial  condition. 


4.4.2  Results. 

Using  the  results  from  STEP  for  a  100  kg  unconfined  explosion  of  TNT  at  a  height  of  four  meters 
as  the  basis  for  this  simulation,  a  buoyant  bubble  of  gas  and  1  kg  of  dust  was  mapped  into 
CATAPULT.  The  simulation  was  conducted  in  axisymmetric  coordinates  and  the  results  were 
mostly  as  would  be  expected.  The  bubble  initially  expanded  slightly  and  then  rose  to  a  height  of 
85  m  while  slowly  dissipating.  All  the  while,  dust  slowly  falls  to  the  ground  beneath  the  initial 
explosion.  Figure  15  shows  a  series  of  plots  of  gas  density  overlain  with  contour  lines  of  the  log  of 
Dust  density.  Given  the  relatively  little  data  available  for  this  type  of  explosion  and  the  generally 
poor  output  available  from  STEP  for  unconfined  explosions,  our  efforts  were  primarily  focused  on 
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(a)  Time  =  7  seconds 


(b)  Time  =  17  seconds  (c)  Time  =  27  seconds 


(d)  Time  =  47  seconds  (e)  Time  =  57  seconds  (f)  Time  =  67  seconds 


(g)  Time  =  87  seconds  (h)  Time  =  97  seconds  (i)  Time  =  107  seconds 

Figure  14.  Gas  density  contours  (shaded)  with  dust  density  contour  lines  over-lain  for  the  3D 
simulation  taken  along  a  y-z  cross-section  from  the  center  of  the  domain. 
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other  phases  of  this  project. 

As  discussed  in  Section  4.2,  several  simulations  of  dust  bubbles  in  axi  symmetric  and  3-D 
coordinates  were  conducted  and  used  as  a  validation  for  the  3-D  component  of  CATAPULT. 
These  results  were  remarkably  similar  to  those  observed  in  the  unconfined  explosion  simulations. 


4.5  Nuclear  Cloud  Rise  Phenomena. 

The  late-time  evolution  of  nuclear  particulate  clouds  were  modeled  at  the  behest  of  the  program 
director.  Initial  conditions  were  provided  from  two  separate  sources:  (1)  Data  from  MAZ 
computations  was  provided  courtesy  of  Titan  Research  Technologies  Corp.  and  (2)  data  from 
SHARC  computations  was  provided  courtesy  of  Maxwell  Technologies,  Inc. 

Two-dimensional  axisymmetric  computations  are  performed  for  a  500  ton  surface  burst  for  a 
computational  domain  of  1000x3000  meters.  The  initial  data  is  mapped  from  a  two-dimensional 
axisymmetric  compressible  computation  from  MAZ.  The  MAZ  computation  simulates  the  cloud 
evolution  for  a  period  of  0-60  seconds  after  burst.  Care  must  be  taken  to  ensure  that  data  used 
from  the  compressible  computation  is  at  an  effectively  incompressible  state  throughout  the  flow 
field.  Three  flow  field  criteria  are  used  to  determine  the  time  when  the  flow  is  sufficiently 
incompressible:  (1)  an  absence  of  shock  waves,  (2)  a  sufficiently  low  Mach  number,  and  (3)  a 
sufficiently  low  velocity  divergence.  Although  the  shock  wave  from  the  initial  blast  leaves  the 
domain  of  interest  at  a  very  early  time,  the  Mach  number  and  velocities  remain  high  for  several 
seconds.  At  20  seconds  after  blast,  the  maximum  Mach  number  present  in  the  flow  field  is  Mmax 
=  0.65  and  the  average  Mach  number  (averaged  for  grid  points  having  M  >  0.1)  is  Mavg  =  0.30. 
At  60  seconds  the  maximum  Mach  number  is  decreased  to  Mmax=  0.48  and  the  average  Mach 
number  is  decreased  to  Mavg—  0.25.  A  contour  plot  of  the  velocity  diverge  calculated  using 
centered  differences  is  shown  in  Figure  16  for  the  20  and  60  second  times.  Figure  16a  indicates  a 
sizable  region  of  compressible  flow  located  in  the  rising  fireball  with  a  maximum  velocity 
divergence  of  V  •  u  =  0. 1635s-1.  At  60  seconds  the  flow  field  is  largely  incompressible  as  indicated 
in  Figure  16b  where  the  maximum  divergence  V  •  u  =  0.035s-1  is  decreased  nearly  an  order  of 
magnitude  from  that  of  Figure  16  a. 

As  expected,  attempting  to  advance  the  solution  with  the  present  method  using  initial  data  that 
has  regions  of  compressible  flow  is  found  to  cause  numerical  instabilities.  The  evolution  of  the 
density  fields  for  both  the  gas-phase  and  particle-phase  computed  from  the  multiphase  projection 
method  are  shown  in  Figure  17  for  the  500  ton  surface  burst  described  above.  The  initial 
condition  is  mapped  from  MAZ  results  at  60  seconds  after  blast  (AB).  Nine  frames  depict  the 
evolution  of  the  nuclear  cloud  for  the  times  indicated  on  each  frame  in  Figure  17.  The  solution  is 
advanced  for  six  minutes,  or  420  seconds  AB.  This  required  approximately  8  hours  of  CPU  time 
on  a  Pentium  II  400  MHz  processor  for  32x128  grid  points.  The  color  scale  depicts  the  density  p 
of  the  gas-phase  which  indicates  the  atmospheric  stratification  of  the  air,  where  redder  shades 
indicate  higher  densities.  The  contour  lines  depict  the  log  of  the  dust-particle  density  p.  The 
specified  particle  diameter  of  60  microns  is  determined  by  calculating  the  mass-weighted  average 
for  eight  particle  groups  present  in  the  upper  cloud  of  the  MAZ  computation.  The  dust  contours 
of  the  initial  condition  depicted  in  the  60  second  frame  of  Figure  17  shows  the  major  structures 
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Figure  15.  Gas  density  contours  (color)  with  dust  density  contour  lines  over-lain  for  the  axisym- 
metric  simulation  of  an  unconfined  explosion. 
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Figure  16.  Velocity  divergence  for  a  nuclear  cloud-rise  simulation  computed  from  MAZ  at  (a) 
20  seconds  after  burst,  and  (b)  60  seconds  after  burst.  The  divergence  at  60  seconds 
has  decreased  by  an  order  of  magnitude  from  that  at  20  seconds  indication  a  flow 
field  that  is  much  more  incompressible. 
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typical  to  nuclear  clouds,  consisting  of  the  pedestal,  stem,  and  the  cloud-top. [1]  The  shade  of  the 
cloud  top  is  much  lighter  than  that  of  the  surrounding  atmosphere,  which  indicates  that  the 
fireball  has  a  lower  density  than  that  of  the  ambient  air. 

The  top  of  the  cloud  initially  rises  due  to  both  the  initial  upward  momentum  of  the  cloud  and  the 
buoyancy  of  the  cloud  created  by  low  density  in  the  cloud-top.  Both  dust  and  the  high-density  air 
at  the  lower  altitude  are  entrained  into  the  toroidal  vortex  as  the  cloud  rises  into  the  stratified 
atmosphere  as  shown  in  the  108  second  frame  in  Figure  17b.  After  reaching  a  maximum  height  of 
around  1900  meters,  the  high-density  air  that  is  convected  to  high  altitudes  by  the  cloud  begins 
fall  back  toward  ground  as  shown  in  Figure  17c.  This  results  in  the  convection  of  low-density  air 
in  the  upper  atmosphere  to  lower  altitudes,  indicated  by  the  gas  density  field  in  Figure  17d.  The 
low-density  air  then  ’’bounces”  back  into  the  upper  altitudes  once  again  entraining  high-density 
air  to  higher  altitudes  and  driving  the  cloud  top  upwards  to  reach  a  maximum  altitude  at  over 
2100  meters  (Figures  17(e-g)).  The  gas  phase  then  begins  a  second  oscillation  as  indicated  in 
Figures  16(h-i). 

While  the  dust-phase  is  shown  to  generally  mirror  the  movement  of  the  gas-phase,  the  dust  is  also 
observed  to  fall  toward  the  lower  boundary  due  the  gravity  force.  This  effect  is  observed  by  the 
movement  of  the  dust  density  contours  near  the  bottom  of  the  stem.  The  dust-phase  in  the 
cloud-top  is  also  observed  to  spread  in  the  radial  direction  via  the  main  toroidal  vortex  as  well  as 
a  secondary  vortex  created  by  the  oscillation  of  the  gas-phase.  Oscillating  phenomena  such  as  that 
observed  in  the  computation  described  above  are  well  known  in  the  field  of  meteorology  as 
atmospheric  gravity  waves.  A  parcel  of  fluid  vertically  displaced  from  its  equilibrium  position  in  a 
region  of  stable  stratification  will  oscillate  at  a  frequency  that  is  a  function  of  the  vertical  density 
gradient.  This  frequency  is  called  the  Brunt-Vasalla  frequency  and  can  be  approximated  as 
follows[23], 
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The  frequency  for  the  present  case  is  calculated  to  be  4.835:cl0~3  Hz,  or  a  period  of  207  seconds, 
using  a  linear  approximation  of  the  density  gradient  over  the  range  of  the  maximum  cloud-height 
and  the  average  density  of  the  gradient  for  Eq.  12.  This  calculated  value  is  in  good  agreement 
with  the  period  of  198  seconds  observed  in  the  computation. 


Sensitivity  of  the  specified  particle  diameter  d  is  demonstrated  by  the  computation  shown  in 
Figure  18.  Here,  the  particle  diameter  for  the  same  500  ton  surface  burst  is  specified  to  be  200 
microns,  determined  by  a  mass  average  of  the  particle  phase  from  the  MAZ  data  that  now 
includes  particles  in  the  pedestal  and  stem  of  the  cloud.  The  200  micron  cloud  in  Figure  18 
evolves  much  in  the  same  manner  as  the  60  micron  cloud  in  Figure  17.  However,  due  to  the 
increased  mass  and  inertia  of  the  individual  particles,  the  cloud  achieves  a  maximum  altitude  of 
approximately  1900  meters,  200  meters  less  than  the  case  in  Figure  2.  Due  to  a  higher  terminal 
velocity,  the  200  micron  particle-phase  is  also  observed  to  fall  out  much  more  quickly  than  the  60 
micron  case  as  indicated  by  absence  of  the  cloud  pedestal  in  Figures  18(d-i)  as  well  as  the 
movement  of  the  dust  contours  in  the  cloud  stem. 


A  nuclear  cloud-rise  simulation  using  SHARC  data  as  the  initial  condition  is  also  performed  to 
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Figure  17.  Nuclear  cloud-rise  simulation  computed  from  the  multiphase  projection  methods  for  £ 
particle  diameter  of  60  microns.  Initial  conditions  are  remapped  from  MAZ  results  at 
60  seconds  after  a  500  ton  surface  burst. 
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Figure  18.  Nuclear  cloud-rise  simulation  computed  from  the  multiphase  projection  method  for  a 
particle  diameter  of  200  microns.  Initial  conditions  are  remapped  from  MAZ  results 
at  60  seconds  after  a  500  ton  surface  burst. 
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(c)  32.6  seconds 


Figure  19.  Nuclear  cloud-rise  simulation  computed  from  the  multiphase  projection  method  for  a 
particle  diameter  of  100  microns.  Initial  conditions  are  remapped  from  SHARC  results 
at  50  seconds  after  a  10  kiloton  surface  burst. 
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demonstrate  the  versatility  of  the  projection  scheme  to  advance  flow  field  solutions  for  a  variety  of 
compressible  data.  Figure  19(a)  shows  the  nuclear  cloud  data  obtained  from  SHARC  for  a  10 
kiloton  yield  at  a  height  of  burst  of  30  meters  23  seconds  after  burst.  The  particle  diameter  is 
specified  to  be  100  microns.  Similar  to  Figures  17  and  18,  the  color  scale  depicts  the  gas  density 
while  the  black  contour  lines  show  the  log  of  the  particle  density.  In  the  same  manner  as  the 
computations  started  from  the  MAZ  data  in  Figures  16,  the  low  density  fire  rises  due  to  both 
buoyancy  and  initial  upward  momentum.  The  cloud-top  exits  the  computational  domain  as  shown 
in  Figure  19(c).  For  this  case,  the  vertical  dimension  of  the  physical  domain  is  too  small  to 
observe  the  oscillation  phenomena  present  in  the  cases  shown  in  Figures  17-18. 

The  case  for  the  500  kiloton  surface  burst  with  60  micron  particles  shown  in  Figure  18  is  mapped 
onto  a  coarse  three-dimensional  grid  and  advanced  out  to  138  seconds  after  blast  shown  in  Figure 
20.  The  initial  flow  field  is  shown  in  Figure  20(a)  where  the  blue  surface  represents  the 
10-06 kg/m3  isosurface  of  the  dust  density  and  the  black  contour  lines  show  the  gas  density.  The 
grid  resolution  is  32x32x64  for  a  domain  of  2000x2000x2500  meters  and  required  approximately  24 
hours  to  run  on  a  400  MHz  Pentium  II  processor.  Similar  to  the  trend  observed  with  the  confined 
plumes  described  in  section  4.3.2,  the  three-dimensional  cloud  reaches  a  maximum  height  much 
lower  than  the  axisymmetric  case,  approximately  1750  meters  shown  in  Figure  20(b)  at  78  seconds 
after  burst  as  compared  to  2100  meters  in  Figure  18.  In  the  same  manner  as  the  axisymmetric 
case  the  high  density  gas  and  entrained  dust  fall  back  toward  the  lower  boundary  shown  in  Figure 
20(c). 
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(c)  32.6  seconds 

Figure  20.  3-D  nuclear  cloud-rise  simulation  computed  from  the  multiphase  projection  method  for 
a  particle  diameter  of  100  microns.  Initial  conditions  are  remapped  from  TRT  data. 
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Section  5 
Conclusions 


The  overall  objective  of  this  effort  was  to  develop,  implement,  and  demonstrate  high  resolution 
computer  codes  to  simulate  multiphase  reactive  flow  fields  associated  with  cloud-rise 
phenomenology.  These  codes  are  important  to  DTR.A  interests  and  are  applicable  to  a  broad 
range  of  DTRA  and  defense  applications,  particularly  post-explosion  environments  for  nuclear, 
biological,  and  chemical  (NBC)  weapons.  Such  codes  also  have  potential  commercial  application 
and  nonmilitary  benefits  since  multiphase  particulate  flows  are  also  encountered  in  numerous 
industrial  settings. 

To  meet  the  overall  objective,  KTI  has  developed  a  high  performance  numerical  modeling  software 
package  named  the  Computational  All-Time  Algorithm  Package  with  Underlying  Lawrence 
Technology  (CATAPULT)  for  the  numerical  simulation  of  compressible  and  incompressible 
multiphase  flows.  ”  All-Time”  indicates  CATAPULT’S  capability  to  model  flows  in  both  the 
compressible  and  incompressible  time  scales.  Having  both  compressible  and  incompressible  flow 
solver  modules  provide  an  efficient  method  to  advance  solutions  for  transient  flows  that  are 
compressible  at  early  times  and  incompressible  at  later  times.  CATAPULT  incorporates  the  latest 
technologies  in  high-order  Godunov  (HOG)  schemes  and  projection  methods  and  has  the 
capability  to  handle  two-dimensional  planar,  axisymmetric,  and  full  three-dimensional 
computations. 

The  development  and  utilization  of  CATAPULT  has  resulted  in  an  increased  understanding  of 
post-explosion  environments  both  in  terms  of  the  physical  processes  involved  and  the  issues  that 
arise  in  developing  accurate  models  for  such  environments.  The  technology  and  knowledge 
resulting  from  this  effort  serve  as  a  significant  step  towards  the  future  development  of 
sophisticated,  accurate,  and  fast  computational  tools  that  will  further  the  knowledge  of 
multiphase, _  particulate  flow  field  environments.  The  following  subsections  give  a  brief  review  of 
the  developments  and  findings  in  this  study,  a  discussion  on  the  extent  to  which  the  objectives 
were  achieved,  and  the  work  required  to  develop  the  technology  further. 

5.1  Results. 


The  CATAPULT  software  was  developed  and  assembled  based  on  the  suite  of  advanced  numerical 
schemes  developed  by  the  research  group  at  LBNL.  The  compressible  flow  module  consists  of  the 
scheme  developed  by  Collins  et  al.  [12]  and  is  tailored  to  simulate  NBC  cloud-rise  scenarios. 

The  incompressible  flow  module  was  developed  that  incorporates  a  particle-phase  model  like  that 
in  Collins  et  al.[12]  with  the  projection  method  of  Bell  et  al.  [5].  The  solution  method  for  the 
governing  gas-phase  equations  is  modified  with  an  intermediate  projection  as  described  in  [3]  and 
the  solution  for  the  governing  particle-phase  equations  is  modified  with  unsplit  Godunov  scheme 
for  the  particle  equations  as  described  in  [9].  The  rate  of  convergence  for  the  resulting  multiphase 
projection  method  indicates  second-order  accuracy  for  the  gas  phase  and  between  first  and 
second-order  accuracy  for  the  particle-phase. 

CATAPULT  was  used  to  investigate  plumes  resulting  from  confined  explosions.  Through 
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comparison  with  experimental  results  and  empirical  models,  it  has  been  shown  that  plume 
phenomena  are  not  well  reproduced  by  axisymmetric  compressible  models.  Even  with  extremely 
high  resolution,  the  basic  physics  of  the  problem  is  not  represented  by  an  axisymmetric 
assumption.  Although  the  axisymmetric  and  3D  results  agreed  well  in  the  near  field  region,  in  the 
far  field  where  the  plume  tip  was  located  after  less  than  a  second,  only  the  3D  simulation  was  able 
to  accurately  reproduce  experimental  results.  The  plume  height  and  spreading  rate  were  both 
captured  in  a  much  better  manner  by  a  fully  three  dimensional  model  at  later  stages  of  the  flow  in 
the  compressible  regime.  The  ensemble  velocity  profiles  of  the  3D  simulation  fit  the  expected 
theoretical  profiles  to  a  reasonable  extent,  but  the  actual  asymmetry  in  the  flow  proved  to  be  a 
major  contributor  to  the  resulting  plume  height. 

The  entrainment  rates  observed  in  the  compressible  phase  of  the  3D  flow  were  slightly  lower  than 
the  rates  reported  by  Bremhorst  and  Hollis  [7]  for  near  field  values,  but  began  to  increase  quickly 
in  the  far  field.  The  linear  trend  of  plume  height  versus  the  square  root  of  non-dimensional  time 
was  observed  during  the  early  phase  of  the  compressible  problem,  but  did  not  persist  until  the 
vent  velocity  reached  zero.  The  findings  of  this  study  indicate  that  the  results  from  previous  jet 
flow  researchers  are  applicable  to  multi-phase  plumes  in  the  near  field  region.  However,  as  the  flow 
expanded  away  from  the  near  field  to  distances  greater  than  100  d0,  the  plume  characteristics 
generally  began  to  diverge  from  jet  characteristics  as  entrainment  rates  increased  and  plume  tip 
penetration  slowed. 

Finally,  the  incompressible  phase  of  the  flow  has  also  proven  to  be  poorly  suited  to  the 
axisymmetric  constraint.  The  axisymmetry  in  the  problem  continues  to  predict  large  scale 
coherent  structures  in  the  flow,  that  are  most  likely  not  present.  In  addition  to  incorrectly 
capturing  the  velocity  fields,  this  tendency  toward  coherent  structures  in  the  axisymmetric  model 
can  significantly  affect  the  distribution  of  dust  in  a  multi-phase  flow.  Plume  height  and  general 
shape  were  also  found  to  be  inadequately  captured  with  an  axisymmetric  model,  especially  when 
run  at  higher  resolutions. 

Multiphase  particulate  clouds  resulting  from  unconfined  explosions  were  also  simulated.  There  is 
very  little  existing  data  for  unconfined  explosions  that  can  be  used  for  comparison  and  as  initial 
conditions  for  the  simulations.  Thus,  the  unconfined  simulations  in  this  effort  are  limited  to  the 
output  of  the  unconfined  explosion  model  component  of  STEP,  which  is  relatively  undeveloped  as 
well.  While  axisymmetric  and  three-dimensional  simulations  appear  to  qualitatively  capture  the 
large-scale  dynamics,  more  detailed  initial  conditions  are  required  to  obtain  more  quantitatively 
accurate  unconfined  simulations. 

Cloud-rises  resulting  from  nuclear  blasts  for  the  incompressible  time-scales  were  simulated  at  the 
behest  of  the  program  director.  Initial  conditions  were  obtained  from  two  different  hydrocodes 
that  have  been  used  to  simulate  the  compressible  time  scale  of  nuclear  blasts.  MAZ  data  was 
provided  courtesy  of  TRT  and  SHARC  data  was  provided  courtesy  of  Maxwell  Technologies.  Data 
from  both  sources  were  successfully  mapped  into  the  incompressible  module  of  CATAPULT  and 
advanced  to  later  (incompressible)  times. 

An  example  of  one  such  case  is  the  computation  of  an  axisymmetric  cloud  evolution  for  a  500  ton 
surface  burst.  The  computation  is  initialized  using  MAZ  flow  field  data  at  60  seconds  after  burst. 
The  60  second  flow  field  is  determined  to  be  incompressible  and,  therefore,  appropriate  to  begin 
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the  projection  calculation.  The  solution  is  advanced  out  to  a  time  of  six  minutes  from  the  initial 
condition  and  accurately  captures  the  large  scale  structure  of  the  cloud  evolution,  including  the 
Brunt- Vasalla  oscillation.  The  maximum  altitude  obtained  by  the  cloud  and  the  rate  of  fall  out 
for  the  particle-phase  is  found  to  be  sensitive  to  the  specified  particle  diameter. 


5.2  Progress  vs.  Proposed  Objectives. 

The  objectives  listed  in  the  Phase  II  proposal  are  restated  below: 

Technical  Objective  1:  Simulate  gas-particle  cloud  rise  resulting  from  the  compromise  of  a 
confined  CB  facility.  Analyze  the  results  and  compare  with  experimental  data  and/or  other 
models  as  available. 

Technical  Objective  2:  Simulate  gas-particle  cloud  rise  resulting  from  the  compromise  of  an 
unconfined  CB  facility.  Analyze  the  results  and  compare  with  experimental  data  and/or  other 
models  as  available. 

Technical  Objective  3:  Modify,  implement,  and  optimize  algorithms,  schemes,  and  codes  and/or 
other  models  as  available. 

The  first  objective  was  achieved  through  the  development  of  CATAPULT  to  simulate  both  the 
compressible  and  incompressible  time-scales  of  a  confined  scenario.  Knowledge  and  understanding 
of  the  physical  processes  and  issues  of  modeling  were  garnered  through  comparisons  with 
experiments,  i.e.,  the  Dipole  East  159  [2]  and  (jet  experiments  from  the  literature),  as  well  as 
results  from  STEP. 

The  second  objective  was  partially  achieved  through  CATAPULT  simulations  for  the 
incompressible  time-scale  of  an  unconfined  event.  Initial  conditions  were  limited  to  the  relatively 
undeveloped  unconfined  explosion  module  of  STEP.  Experimental  data  was  not  available  for 
comparison. 

The  third  objective  was  fully  achieved  through  the  developments  of  CATAPULT  described  in  the 
first  two  paragraphs  in  the  previous  sub-section. 


5.3  Future  Work. 

The  CATAPULT  code  developed  under  this  contract  is  a  significant  addition  to  the  modeling 
capabilities  of  DTRA.  In  order  to  become  a  useful  operational  tool  however,  some  further  work  is 
required.  The  dust  phase  of  flow  in  CATAPULT  is  currently  assumed  to  consist  of  dust  particles 
of  constant  size.  Modeling  the  dust  using  multiple  particle  sizes  would  increase  the  capabilities  of 
the  code  significantly.  This  could  be  accomplished  in  several  ways  and  some  study  would  be 
required  to  determine  the  most  effective  method.  Also,  in  order  to  correctly  model  nuclear  cloud 
rise  phenomena,  the  addition  of  thermodynamic  coupling  to  the  code  would  be  required.  Most 
nuclear  clouds  contain  large  amounts  of  water  vapor  that  cool  and  condense  as  it  is  entrained  to 
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higher  altitudes.  This  condensation  releases  a  significant  amount  of  latent  heat  that  contributes 
further  to  the  energy  and  buoyancy  of  the  cloud.  In  order  to  model  nuclear  cloud  rise  events  in 
most  environment,  it  is  essential  that  thermodynamics  be  included  in  the  model.  Other  physics 
prominent  in  NBC  post-explosion  environments  should  also  be  incorporated  into  the  model 
including  non-dilute  regions  of  particulates  and  chemical  reactions. 

Perhaps  the  most  important  improvement  required  to  make  CATAPULT  a  useful  operational  tool 
is  code  parallelization.  Three  dimensional  computations  of  plume  phenomenon  from  confined 
explosions  required  approximately  two  months  to  run  to  completion  on  Cray  vector  computers 
and  PC’s.  Code  parallelization  would  allow  CATAPULT  to  run  a  similar  simulation  to  completion 
in  12  to  36  hours  on  a  typical  Cray  T3E  machine,  and  make  operational  use  of  the  code  highly 
feasible.  The  technology  necessary  to  accomplish  code  parallelization  currently  exists  and  is 
relatively  uncomplicated  to  implement,  although  it  does  require  a  significant  effort  in  terms  of 
manpower  and  computer  resources  to  test  and  validate. 
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